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ABSTRACT 

Accretion discs at sub-pc distances around supermassive black holes are likely to cool rapidly 
enough that self-gravity results in fragmentation. Here, we use high-resolution hydrodynamic 
simulations of a simplified disc model to study how the outcome of fragmentation depends 
upon numerical resolution and cooling time, and to investigate the incidence of binary for- 
mation within fragmenting discs. We investigate a range of cooling times, from the relatively 
long cooling time-scales that are marginally unstable to fragmentation down to highly unsta- 
ble cooling on a time-scale that is shorter than the local dynamical time. The characteristic 
mass of fragments decreases with reduced cooling time, though the effect is modest and de- 
pendent upon details of how rapidly bound clumps radiate. We observe a high incidence of 
capture binaries, though we are unable to determine their final orbits or probability of sur- 
vival. The results suggest that faster cooling in the parent disc results in an increased binary 
fraction, and that a high primordial binary fraction may result from disc fragmentation. We 
discuss our results in terms of the young massive stars close to the Galactic Centre, and sug- 
gest that observations of some stellar binaries close to the Galactic Centre remain consistent 
with formation in a fragmenting accretion disc. 
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1 INTRODUCTION 

In recent years the topic of gravitational instabilities in accretion 
discs has been the subject of a great deal of theoretical and nu- 
merical research. Such gravitational instabilities have applications 
both to accretion disc theory, in terms of the transport of angular 
momentum, and also to theories of star and planet formation, as 
under ce rtain circumstances the instability can lead to fragmenta- 
tion (e.g. lDurisen et al.ll2007l : lLodato 2007) . In particular, the frag- 
mentation of accretion discs has b een invoked to explain both th e 
formation of gas giant planets (e.g. lBosil 19971 : iMaver et alj|2002h . 
and also to explain the formation of young, massive stars observed 
close to super-massive black holes at the centres of galaxies (e.g., 
IPolvachenko & Shukmanll 19771 : iLevin & Beloborodov|[2003l) . 

In order for a gravitationally unstable disc to fragment, two 
conditions must be met. The first is that the disc must be sufficiently 
cold, or sufficiently massive, to satisfy thc iToomra ( il964) criterion 



Q 



TlGI. 



< 1 



(1) 



Here is the local sound speed, k the epicyclic frequency (which 
in a thin disc is approximately equal to the local orbital frequency), 
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and E the disc surface density. In addition, in order to fragment 
the disc must cool suffi ciently rapidly to overcome compressional 
heating during collapse. iGammiel ( 1200 ih used a local analysis, and 
found the fragmentation boundary to be 

'cool < yScnti^"' , (2) 

where f^ooi is the local cooling time-scale, Q is the orbital fre- 
quency, and jScrii - 3. Subsequent studies u sing global simulations 
have verified the validity of this crit erion jRice et al]|2003l , [20051 : 
iMeifa et alj|2005l : iBolev et alj|2007l) . although the exact value of 
jScrit required for fragmentation can be larger by a factor of 2- 
3 depending on the details of the thermal physics adopted (cool- 
ing law, equation of state, etc.). To date, however, most work 
has either made use of simple th ermal physics (e.g. iRice et aD 
I2OOI I2OO5I : lNavakshinet"ai]|2007h . or used more realistic thermal 
phys ics in the specific case of self-gravitatin g protop lanetary discs 
(e.g. lCai et al.ll2007l : [st"amatellos & Whitworthll2008h . More gener- 
ally, most studies of fragmentation have tended to look at cases of 
critical cooling (with Clt^oo] - Pent), despite the fact that some discs, 
notably thos e around black holes, are thought to cool much more 
rapidly (e.g. lBertin & Lodato 2001; Goodman 2003). 

In this paper we use numerical hydrodynamics to study the 
influence of rapid cooling on the fragmentation properties of grav- 
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itationally unstable accretion discs. In Section|2]we present thie de- 
tails of our numerical simulations, in which we adopt a cooling law 
consistent with radiative cooling in an optically thick disc. We use 
this model to study the non-linear evolution of the instability over 
a wide range in cooling times. We also consider in detail the effects 
of numerical resolution on our calculations, varying the number of 
fluid elements by a factor of 64 between our lowest- and highest- 
resolution calculations. We present our results in Section [3] first 
considering issues of numerical convergence and then addressing 
the effects of varying the cooling time. In Section|4]we discuss the 
limitations of our analysis, and the implications of our results for 
the formation of stars near to the Galactic Centre, before summa- 
rizing our conclusions in Section[5] 



2 NUMERICAL METHOD 

The simulations presented in this paper use the publicly-available 
smoo thed-particle hydrodynamics (SPH) code gadget2 (Springel 
l2005h . We have modified the code to include a simple cool- 
ing prescription (discussed below), which is representative of the 
cooling expected in a real ac cretion disc. We adopt the standard 
iMonaghan & Gingoldl l ll983h prescription for t he SPH visco sity 
(with QfspH = 1-0), using the "Balsara-switch" l lBalsaralll995h to 
limit the artificial sh ear viscosity (as specified in Equations 1 1- 
12 o f lSpringelll200^ . We allow for a variable gravitational soften- 
ing length, and fix the SPH smoothing and gravitational softening 
lengths to be equal throughout. We use the standard Barnes-Hut 
formalism to c ompute the gravitational force tree (as described in 
ISpringell200^ . and use A'ngb = 64 ± 2 as the number of SPH neigh- 
bours. Typically, a handful of SPH particles end up in rather close 
orbits around the central gravitating mass and, if left unchecked, re- 
sult in unreasonably short time-steps. Consequently we use a single 
sink particle for the central gravitating mass. The sink particle sim- 
pl y accretes all gas pa rticles that pass within its radius (as described 
m ICuadra et al.ll2006 ll. and we set the sink radius to be 1/4 of the 



inner disk radius. This is simply a numerical convenience adopted 
in order to save computing time, and has no physical effect on the 
simulations. The simulations were run on the tungsten Xeon Linux 
cluster at NCs41 

using 64 parallel CPUs for the highest-resolution 

runs. 

Numerical simulations of gas disc dynamics suffer from a ba- 
sic computational limitation, namely that the pressure scale-height 
of the disc acts as a minimum length-scale that must always be 
well-resolved. Requirements on computing time are therefore very 
strongly dependent on the disc aspect ratio (the required time typ- 
ically scales as (H/Ry^ or steeper), and accurate global simula- 
tions of thin discs are very computationally expensive. In a self- 
gravitating disc with Q - I, the disc aspect ratio is approximately 
equal to the ratio of the disc mass to the central object mass, and 

relatively 

_ t alj|2003l 

12001 iBolev etai]|2007h . This is typical of young protoplanetary 
discs, but discs around black holes may in fact be much thinner 
(H/R < 10"^). We therefore seek a method of simulating the be- 
haviour expected in a very thin disc, but without incurring the se- 
vere computational penalty that usually results. 

We make use of an artificial disc model, which behaves like a 
thin self-gravitating disc while in fact remaining thick enough to be 



' See |http : //www . ncsa . uiuc . edu/| 



to date most numerical simulations have focus sed on r elatively 
massive, thick discs (w ith Mj^c ~ 0. IM.. e.g., iRice etl 



well-resolved even at moderate resolution. Spiral density waves in 
self-gravitating discs tend to be dominated by a relatively narrow 
range in spiral mode number, with progressively thinner discs re- 
sultin g in ever higher-order spiral waves (see the review bv lLodatol 
|2007). Moreover, as the disc thickness is reduced, the "local ap- 
proximation" (that global modes of the instability are negligible) 
becomes more accurate. Consequently, we adopt a relatively large 
disc mass, but impose an artificial radial cooling profile that results 
in the disc fragmenting in a single, near-circular filament. The anal- 
ysis of our simulations considers only this filament, with the rest of 
the disc regarded as a boundary condition: the dynamics of frag- 
mentation in a single, near-circular filament are a good approxima- 
tion to a real self-gravitating disc in the thin disc limit. With this set- 
up we are able to resolve the disc scale-height properly with only 
moderate numbers of SPH particles (see Section [23] below), and 
can thus run simulations with much larger particle numbers ("over- 
resolving" the disc thickness) to study the fragmentation process in 
detail when the cooling time is short. 



2.1 Initial Conditions 

We adopt initial conditions similar to those that h ave been adopted 
in previous studies of disc fragmentation (e.g. iRice et al] l2003l 
I2OO5I : [Alexander et ar]|2008D . We adopt a surface density profile 



E(r) oc 

and a disk temperature that scales as 

T(r) oc r^l^ , 



(3) 



(4) 



where r is the (cylindrical) radius. The disc sound speed c, oc r''^ 
is normalized so that the Toomre parameter (Equation[T) is equal to 
2.0 at the outer boundary. The disc is initially vertically isothermal, 
so the vertical structure of the disc is given by 



S(r) 



exp 



z 



(5) 



where H = cJCl is the disc scale-height (here Q = ■\/GMJr^ is 
the Keplerian orbital frequency.). With this set-up the initial disc is 
marginally gravitationally stable, and the instability is allowed to 
develop in a physical manner as the disc cools. We define a disc 
mass of Md = 0. IM, (where M, is the mass of the central gravi- 
tating mass), and model the disc using A'sph SPH particles of mass 
m = Mi/Nsm- We adopt a system of units where M, = 1, the inner 
edge of the disc is at r = 1, and the time unit is the orbital period at 
r = 1. (Consequently, G = An- in code units.) The initial distribu- 
tion of the SPH particles is obtained by randomly sampling the disc 
density profile described above, with each particle given a circular 
velocity in the x-y plane. We adopt an adiabatic equation of state 
throughout, with adiabatic index 7 = 7/5 (as expected for a disc 
consisting primarily of molecular hydrogen). 

2.2 Cooling 

The principle aim of this study is to investigate the effect of variable 
cooling rates on the fragmentation properties. At each time-step, we 
allow the internal energy of the ith particle, u,, to cool as 



ldu\ Hi 

\ /cool., 'mol 



(6) 



The cooling time-scale tcao\ is defined in terms of the dimensionless 
parameter /? and the orbital frequency 
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(7) 



As mentioned in Section [T] previous studies have found that the 
fragmentation boundary lies in the range /3crit - 3- 7.5, with a weak 
depe ndence on the adopted equation of state (e.g. Gam mie|l200ll : 
iRice et al. 2003 . .2005) . Here, in order to study the effects of vari- 
able cooling, we define yS to be a function of radius, density and 
time: 



P = 'R(r)D(p)T(t). 



(8) 



As discussed above, the radial dependence of /? is chosen to mimic 
the behaviour of ajightly-wound spiral density pattern (as expected 
in a thin disc; e.g. lLodatol2007D . but without the excessive demands 
on computing resources that usually arise when simulating very 
thin discs. We choose 7? so that the disc fragments in a ring at r = 2, 
and can therefore save computing power by considering a limited 
range in radius. The only requirement on the radial range is that the 
boundaries of the disc do not influence the fragmentation, so we fix 
the inner boundary at r = 1 and the outer boundary at r = 3. The 
radial dependence of the cooling time is chosen to be a Gaussian 
profile centred on ro = 2 



Kir) = 8.0 X 



, 1 / (r-rof 



(9) 



with A = 0.25. Thus Kir) varies from 4-8, with the minimum value 
at r = ro. The normalisation parameters are chosen so that the disc 
is initially marginally stable to fragmentation everywhere except 
close to r = ro, and the value of A is chosen so that the radial length- 
scale over which the cooling time varies is a few times larger than 
the pressure scale-height of the disc. 

In addition, we allow the cooling time to vary with the local 
density. In a real self-gravitating disc we expect radiative cooling 
to be optically thick (see the discussion in Section|4j, so the cool- 
ing rate is density-dependent. If cooling is optically thick then the 
radiative cooling time-scale is proportional to the photon diffusion 
time-scale. In the case of constant opacity k, the time, to, for a pho- 
ton to diffuse a length L is given by 



to = PK, 

C 



(10) 



where c is the speed of light. In a collapsing (thinning) disc the 
appropriate length-scale is the disc thickness h, and the density p ^ 
I./h. Thus, for 2-D collapse we have 



(11) 



and we see that the cooling time-scale decreases with increasing 
density p. By contrast, in a collapsing spherical clump of mass in 
and radius h, we have p m/h^, so 



to oc h 



(12) 



and therefore the cooling time-scale increases (i.e. cooling becomes 
less efficient) as the clump becomes more dense. This increase in 
the cooling time-scale can lead to significant suppression of cooling 
as the forming fragments reach high density, and consequently we 
adopt a density-dependence with the following form: 



D(p) 



1, - 



1/3 



(13) 



Thus for densities below po the cooling time is independent of the 
density; at higher densities the cooling time increases as p''^. Em- 
pirically, we choose po = I.O in code units (central object mass / 
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Table 1. List of .simulations run, showing resolution and physical param- 
eters for each simulation, ffmj denotes the time (in code units) when the 
first fragment forms, and /?'(?fiag) indicates the approximate value of the lo- 
cal cooling pai'ameter /}'(tfrng) = ^(''o)7'('frag) (i-e- neglecting the density 
dependence) when fragmentation occurs. 



inner disc radius'), as this is typical of the local density in our sim- 
ulations at the point where bound, roughly spherical objects begin 
to form. In practice the weak dependence of the cooling time on 
density means that our results do not depend significantly on the 
exact value of po. 

Lastly, in order to study the effect of very rapid cooling we 
allow the cooling rate to vary with time in the simulations. The 
motivation for this is that one cannot simply begin a numerical sim- 
ulation with a cooling time that is much shorter than the dynami- 
cal time-scale, as the resulting collapse and fragmentation merely 
amplifies noise in the initial conditions and has little or no physi- 
cal meaning. Instead, one must allow the instability to develop in 
a physical manner before the cool ing time is allowed to become 
short. We follow the approach of Iciarke et alj ( |2007|) , and allow 
the cooling time to decrease as 



T(t) = 1 - I - 



(14) 



Here Vb is the local orbital period, and t is a dimensionless pa- 
rameter that we use to parametrize the cooling rate. We note that 
the form of T(t) allows the cooling time to become very small (and 
indeed negative) at large t; in practice, however, our simulation run 
for short enough times that this is not a problem. We adopt values of 
the cooling parameter t = 2, 3, 5, oo (in the latter case the cooling 
time is constant with time). With this set-up the disc is marginally 
unstable to fragmentation at r = 2 for t = oo, while smaller val- 
ues of T result in cooling times that are up to order of magnitude 
shorter. At the end of our simulations the smallest values of p are 
in the range 0.2-0.3 (see also Table[T). 



2.3 Simulations 

The issue of numerical resolution is critical when considering cal- 
culations of fragmentation, as insufficient resolution can result in 
spurious suppression or amplification of fragmentation. We have 
conducted a number of calculations using a range of numbers of 
SPH particles for different values of t (see Table[T]l. Resolution re- 
quirements fo r SPH simulations of fr agmentation were presented 
and tested by ' Bate & BurkertI ( 1 1997b, and extended for the spe- 
cific case of self-gravitating discs by iNelson (2006). In this case 
the Jeans mass must be resolved into > lOA'ngb particles, and the 
disc scale-height H must be resolved into at least 4 SPH smooth- 
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ing lengths in the saturated state of the instabihtjU. For most of the 
models considered here the second of these is the more stringent 
condition, but correctly resolving the Jeans mass becomes more de- 
manding in the simulations with r = 2 (see Section [3Tt . We choose 
our lowest resolution so that we marginally satisfy these conditions 
for the case t = oo, and thus choose Wsph = 200, 000. We then 
move to higher resolution, first increasing the particle number by a 
factor of 8 (to A'sph = 1 .6x lO*"), and then by a further factor of 8 (to 
^SPH = 1.28x 10'). In aperfectly spherical configuration increasing 
^SPH by a factor of 8 results in a factor of 2 (= 8''') increase in the 
spatial resolution. However, in the disc geometry considered here 
we find that a factor of 8 increase in A^sph in fact improves the spa- 
tial resolution by a factor of approximately 2.5. For convenience, 
we subsequently refer to the simulations with A'sph = 200, 000 as 
"low resolution" (LR), A'sph = 1.6 x 10* as "medium resolution" 
(MR), and A'sph = 1.28 x 10' as "high resolution" (HR). For rea- 
sons of computation time we ran the HR simulations only for the 
extreme cases t = oo and t = 2; the LR and MR runs simulations 
were conducted across the full range in r. 

An additional benefit of moving to resolutions higher than the 
minimum requirement is the potential to resolve binary and multi- 
ple systems. Roughly speaking, we expect that individual objects 
will only be able to fragment into multiple systems if they cool 
on a time-scale shorter than the local dynamical time-scale. In a 
self-gravitating disc with Q — i the dynamical time-scale of a sin- 
gle fragment is, to first order, equal to the local orbital time-scale. 
Consequently, we do not expect the formation of "fission" binaries 
in discs with y6 > 1 (although binaries may still form by dynami- 
cal capture events), as thermal pressure acts to smooth out density 
fluctuations on a time-scale shorter than the cooling time-scale. 
Thus the minimum resolution requirement is sufficient in simula- 
tions with p > 1. However, with faster cooling it may be possi- 
ble to form binary and multiple systems via the dynamical fission 
of single clumps. Consequently our high-resolution simulations act 
not only as convergence tests, but also provide a means of studying 
the small-scale fragmentation that may occur in discs with short 
cooling times. 

In order to compare the fragmentation in different simulations 
we must first identify individual fragments ("clumps") in our sim- 
ulations, and we do this simply by adopting a density threshold 
of p = 10.0 (i.e. ten times larger than the density above which 
the cooling time begins to slow). We identify clumps as objects 
which have peak densities above this threshold, and estimate clump 
masses by moving outward in spherical shells one (minimum) SPH 
smoothing length in radius, until the mean density drops below the 
threshold value: all particles within this radius are then deemed to 
be part of the clump. 

One consequence of using a fixed density threshold in this 
manner is that clumps are identified at somewhat earlier times in 
simulations with higher resolution, as the central regions of individ- 
ual clumps are no longer "smeared out" over as long a length-scale. 
This is an artefact of our analysis, rather than a defect in the simula- 
tions, and we solve this simply by comparing simulations at differ- 
ent resolutions relative to the time at which the first clump is identi- 
fied. In practice this difference is very short: the delay in / between 
the identification of the first fragment at the highest and lowest res- 



^ Note that the definition of the SPH smoothing length in lSpringej j2005l) 
differs from that in most of the SPH literature, being larger by a factor 
of 2. Throughout this paper we discuss only the standard definition of the 
smoothing length (e.g. lMonaghaj[l992h . 



Figure 1. Plots of the disc surface density in the simulations with the slow- 
est (t = oo; top row) and fastest (t = 2; bottom row) cooling times, using 
different numerical resolutions. From left-to-right, the three columns show 
the LR, MR and HR simulations respectively, and each plot is made 1/8 of a 
local orbital period (0.35 time units) after the first fragment forms. For slow 
cooling, just below the fragmentation threshold, we see that the fragmenta- 
tion process is well-resolved in all our simulations, although more substruc- 
ture is apparent in the higher-resolution simulations. With fast cooling the 
disc is colder when it fragments, and higher resolution is needed in order 
to resolve the disc scale-height: in this case we see that the LR simulation 
is somewhat under-resolved (e.g. the density peaks at r =: 1.8 are numer- 
ical artefacts), but that A'sph = 1.6 x 10* appears sufficient to resolve the 
fragmentation process properly. 



olutions is a few percent of the total simulation time (see Table[Tll. 
Overall this is a rather crude method for measuring the clump mass 
function, and a more sophisticated approach could, for example, 
allow for non-spherical clumps, or test every particle in order to 
determine how many individual parti cles are bound to each den- 
sity peak (e.g. [Alexander et"ai]|2008l) . In practice, however, such 
"improvements" create at least as many problems as they solve, be- 
cause the treatment of, for example, binary and multiple systems, or 
"circumstellar" discs, becomes ambiguous. The method used here 
is robust across different resolutions, and is sufficiently accurate 
for our needs: a posteriori checks show that negligible fraction of 
clumps identified in this manner are not in fact bound objects 



3 RESULTS 

3.1 Numerical convergence 

We first seek to compare our results at difference numerical reso- 
lutions, in order to test convergence in our simulations. The non- 
linear development of gravitational instabilities essentially ampli- 
fies seed fluctuations in the density field, and it is not possible to 
construct identical noise fields at different particle numbers. Con- 
sequently we instead seek statistical convergence, measuring when 
the first bound objects appear, how many such objects we find, and 
how much mass is present in these objects. We identify clumps us- 
ing the method described in Section [23] and use the mass functions 
derived in this manner to compare different simulations. 

The results of our numerical convergence tests are shows in 
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t (after first frogment) / Inner orbital period 

Figure 2. Results of numerical convergence tests. The upper panel shows 
the evolution of the number of clumps for the simulations with r = oo. The 
origin of the time axis is the time when the first clump is identified, and 
the LR, MR and HR simulations are denoted by black, red and green lines 
respectively. Although more binaries are resolved at higher resolution the 
total mass in clumps does not vary significantly between the runs, indicating 
a good degree of numerical convergence. The lower panel shows the same 
plot for the simulations with t = 2. In this case the LR simulation is rather 
under-resolved, but good convergence is found between the MR and HR 
simulations. 

Figs[I]and|2l In the case of slow cooling (t = oo) we see that 
the fragmentation process is well-res olved in all of our sim ula- 
tions. In terms of the criteria laid out bv lSate & BurkertI ( 1 19971) and 
iNelsonI (l2006h . we resolve the disc scale-height at the mid-plane 
into approximately 5.3 SPH smoothing lengths in the LR run, 12.5 
smoothing lengths in the MR run, and 29.6 smoothing lengths in 
the HR run. In addition the characteristic fragment mass is well- 
resolved throughout so, as noted in Section 123] the LR simula- 
tion marginally satisfies the minimum resolution requirement. Runs 
with higher resolution reveal more substructure, resolving binaries 
and "circumstellar" discs, but fragmentation is independent of res- 
olution. Subsequent interactions between clumps, however, are not 

^ Visualisat ions of our SPH simulations were created using splash: see 
IPricd j2007l) for details. 



independent of resolution: this is clearly seen in Fig[2l where the 
number of clumps beyond At ^ 0.3 differs significantly between 
the LR and HR runs. This is due to mergers of binary clumps, and 
indicates that the interactions between such binaries are not prop- 
erly resolved. This occurs because although the global disc is well- 
resolved the many "circumstellar" discs are not, so the time-scale 
on which binary pairs merge is governed by the rate of numerical 
dissipation of energy in these "sub-discs". Consequently we do not 
attach great significance to our results beyond the time at which 
binaries start to merge. Before clump mergers start to become sig- 
nificant the total mass in clumps varies by less than 20% between 
the LR and HR runs, indicating that good numerical convergence 
has been achieved. 

In the case of fast cooling (t = 2) the disc is colder when it 
fragments, so higher resolution is needed in order to resolve the 
Jeans mass. Here, the pressure scale-height of the disc is again 
marginally resolved in the LRrun (into 4.5 SPH smoothing lengths 
at the disc mid-plane), but the characteristic fragment mass is 
- 4N„gj,m. Thus in this case, as expected, we see that the LR sim- 
ulation is slightly under-resolved (by a factor of approximately 2 
in particle number), but both resolution requirements are satisfied 
comfortably in the MR and HR runs. Binary mergers are somewhat 
more significant in the MR simulation than in the HR simulation 
(again due to under-resolving "circumstellar" discs), but otherwise 
there is good agreement between the MR and HR runs. In both the 
T = oo and T = 2 cases good statistical convergence is found be- 
tween the MR and HR runs, so we conclude that A'sph = 1 .6x 10* is 
sufficient to resolve the fragmentation process properly throughout 
all our simulations. 

3.2 Cooling 

In Section [37T| we found that the MR runs were sufficient to resolve 
the fragmentation process properly across the full range of cooling 
times. We now compare the MR runs across the full range of cool- 
ing times, in order to study the effect of the cooling rate on the frag- 
mentation process. The results of the MR runs are shown in Figs[3}- 
El Fig[3] shows the evolution of the disc surface density in the four 
different simulations, and it is clear that the cooling rate has a dra- 
matic effect on the fragmentation process. Slow cooling leads to 
the formation of larger, more massive fragments with larger separa- 
tions, while faster cooling results in a larger number of less massive 
fragments, with smaller separations. Faster cooling also leads to the 
formation of more binary and multiple systems, although numeri- 
cal dissipation results in many of these systems merging rapidly 
in our simulations. Additionally, the fragmentation process appears 
to be more dynamic in the case of slow cooling, with individual 
objects showing greater perturbations from circular orbits. This oc- 
curs because in the case of slow cooling the collapsing filament 
has time to become unstable to bending modes before it fragments, 
while with faster cooling the filament fragments before it becomes 
significantly bent. 

This qualitative assessment is seen more quantitatively in 
Fig|4l which shows the number and total mass of the bound frag- 
ments in each simulation. Although mass is accreted into fragments 
more rapidly when the cooling is fast, the total mass in fragments 
does not differ significantly with different cooling times. However, 
the simulations with faster cooling produce many more fragments, 
with correspondingly lower masses, and the incidence of multiple 
systems is greatly increased. 

Fig[5] shows how the clump mass function (MF) varies as a 
function of the local cooling time. We compare the MFs at a point 
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Figure 3. Plots of the disc surface density in tlie MR simulations, showing how fragmentation develops as a function of cooling time. From left-to-right, the 
four columns show the evolution of the simulations with r = oo, 5, 3 and 2 respectively. The plots are made in a reference frame that is co-rotating with the disc 
at r = 2, so that only relative motions are evident in the time evolution. (In this frame the disc shear flow appears to move upwards for x <2 and downwards 
for X > 2.) It is clear from the snapshots that faster cooling results in the formation of more, smaller fragments, as well as a higher incidence of multiple 
systems. Note also the large number of binary mergers seen between the second and third snapshots in the simulations with fast cooling. 



where an approximately equal fraction of the disc mass is bound 
into clumps in each simulation, and before the (under-resolved) 
binary mergers begin to dominate the simulations with the fastest 
cooling. We therefore compare the MFs at 0.7 local cooling times 
after the formation of the first fragment (see Fig|4l, although we 
note that the MF comparison is qualitatively similar throughout the 
simulations. All of the MFs are reasonably well-fit by log-normal 
functions, in agreement with previous studies of disc fragmenta- 
tion (e.g. Nayakshin et al. 2007; Alexander et al. 2008). We find 
that faster cooling results in a smaller characteristi c fragment mass. 
This i s in agreement with the previous study of iNavakshin et all 
( l2007l) . which also found that the characteristic fragment mass de- 
creased as the cooling rate increased. In our simulations the local 
(density-independent) cooling time when fragments form varies by 



a factor of ^ 7 between the slowest- (t = oo) and fastest-cooling 
(t = 2) models, and the characteristic mass of the respective MFs 
differs by a factor of 4. We also find that the MFs become some- 
what broader with shorter cooling times, but given the relatively 
small numbers of clumps formed in the t = oo and t = 5 runs 
this result is not statistically significant. The range in characteristic 
masses we find is much smaller than that seen by INavakshin et al] 
( l2007l) . who found that the characteristic mass differed by a factor 
of 100 between their simulations with /? = 3 and /? = 0.3. We 
attribute this to the diff erences between t he coo ling laws adopted 
here and those used by INavakshin et al.l j2007h . In particular, we 
note that the density dependence of our cooling law (which sup- 
presses cooling at high densities) was not included in the simu- 
lations of iNavakshin et al., t2007,) . and the likely consequence of 
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this omission is to enhance fragmentation as the density in clumps 
increases. We therefore suggest that the variation in characteristic 
mass wit h cooling time in real d iscs is less pronounced than that 
found by iNavakshin et al. and conclude that large varia- 

tions in the cooling time-scale result in only modest variations in 
the clump MR 

As mentioned above, a further difference between the fast and 
slow cooling simulations is in the fraction of multiple and binary 
fragments that are formed. The binary fraction is difficult to quan- 
tify due to the rapid merging of the binaries that form, but the dif- 
ferences are clearly seen in Fig[3] With critical cooling (t = oo) 
most objects form as singles, and although some binaries form due 
to dynamical capture events the overall binary fraction is low. How- 
ever, when the cooling is faster the disc is thinner (colder) when it 
fragments, and consequently the clumps are more closely-spaced 
when they first form. This leads to many more capture events, and 
we see many more binaries and multiples in the simulations with 
faster cooling. Indeed, in the simulations with r = 3 and t = 2 
almost all the clumps appear in binary or multiple systems. We see 
no fragmentation of single clumps after the initial fragmentation 
occurs, and all of the binary and multiple systems observed form 
via dynamical capture events (in which the disc gas acts as a "third 
body"). We attribute the absence of any "fission" binaries to the 
density dependence of the cooling time, as the increase in p with in- 
creasing clump density means that the value of fi in the high-density 
clumps rarely, if ever, falls below unity. In this case (as explained 
in Section |23t we do not expect to see sub-clump fragmentation, 
and it seems that the formation of binaries by fission is unlikely in 
real discs with optically thick cooling. The high incidence of cap- 
ture binaries can easily be understood, as faster cooling leads to 
a lower temperature, and therefore a smaller disc scale-height H, 
in the fragmenting disc. The length-scale that is most unstable to 
gravitational instabilities is =^ H, so the typical fragment mass is 
n'LH^ and the typical fragment separation is ^ H. Faster cool- 
ing results in a smaller value of H, which in turn results in smaller 
fragment masses and smaller separations between fragments. In our 
simulations these binaries tend to merge rapidly due to the high rate 
of numerical dissipation in their "circumstellar" discs (as seen be- 
tween f = 5.1 and t = 5.5 in the t = 2 simulation in Fig|3), but 
in reality such mergers would likely be much less efficient. Con- 
sequently, we find that fast cooling in fragmenting self-gravitating 
discs leads to a slight decrease in the characteristic fragment mass, 
and a dramatic increase in the incidence of binary and multiple ob- 
jects. 
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Figure 4. Fragmentation as a function of cooling time. The upper panel 
shows the number of clumps as a function of time in each of the four MR 
simulations, while the lower panel shows the total mass in the clumps. The 
origin of the time axis is set at the point where the first fragment forms, 
and time is plotted in units of the local cooling time (i.e. y3'(ffrag)/f2(ro), 
with the values of /3'(?frag) taken from Table[T). The total mass in clumps is 
independent of the cooling time, but we see that faster cooling results in the 
formation of more objects (with lower masses). 



4 DISCUSSION 

Our simulations are robust and well-understood from a numerical 
perspective, but concerns remain when we consider the physical 
interpretation of our results. The most obviously unphysical aspect 
of our simulations is our "thin disc approximation", using an ar- 
tificial cooling law to allow fragmentation characteristic of a thin 
disc to be modelled at moderate numerical resolution. It would be 
preferable to construct global models of thin discs, but this is not 
computationally feasible with current techniques, especially in the 
case of short cooling times. Our disc model is effectively a local 
approximation, so the primary concern regards the boundary con- 
ditions. In our models the real physical boundaries (at r = 1 and 
r = 3) are sufficiently far from the region of interest that they have 
no influence on the results, and the radial variation of the cooling 
time is negligible in the fragmentation region. However, in a real 



thin disc there are many such filaments, and the typical radial sepa- 
ration between filaments can be as little as a few times H. Dynam- 
ical interactions between adjacent filaments seem likely, especially 
once bound, point-mass-like objects have formed, and the likely 
outcome of this increased rate of clump-clump interactions is the 
formation of even more capture binaries than seen in our simula- 
tions. In some of our simulations (those with slower cooling) some 
of the bound clumps undergo significant radial scattering even in 
the absence of adjacent filaments (moving as far as Ar ^ 0.3 from 
their initial location), and consequently their cooling times change 
in an unphysical way. However, this only occurs at late times in 
the simulations (to which we do not attach great significance), and 
moreover in such dense clumps the dominant term in the cool- 
ing time is the density dependence. Consequently we are satisfied 
that our artificial cooling model is robust, and although our sim- 
ulations should be regarded as numerical experiments rather than 
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Mass function vorlation with cooling rate 




Figure 5. Mass functions of clumps in simulations with different cooling 
times. The mass functions were computed 0.7 local cooling times after the 
formation of the first fragment (i.e. at 0.7 on the time axis in FigE). The 
black/grey, red, green and blue MFs denote the simulations with r = oo, 
5, 3 and 2 respectively. The dotted line denotes the approximate (mass) 
resolution Umit of \ON^^\,m\ all of our clumps are above this limit by at 
least a factor of 5-10. We see clearly that faster cooling results in a smaller 
characteristic mass, and a somewhat broader mass function. 



physical simulations, we are confident attaching limited physical 
significance to our results. 

An additional concern regards the density dependence of our 
cooling law. Fixing the cooling time to scale as p''^ is correct only 
if the opacity is constant and the disc is optically thick, but it is not 
clear that these assumptions hold in reality. In real self-gravitating 
discs the temperatures are expe cted to be ~ lOOK, so dust grains 
domi nate the disc opacity (e.g. iBell & Lin|[l994l : ISemenov et al.l 
l2003h . In the case of a fragmenting disc with Q — 1, the dis c 
is optically thick if the orbital period is < lOOOyr (lLevinll2007h . 
Real discs are generally only thought to be self-gravitating in their 
outer regions, at radii > 50AU i n the case of protostellar discs (e.g. 
IStamatellos & Whitworthll2008l) or at r adii > 0.1 pc in the case of a 
disc around a ~ IO^Mq black hole (e.g. lKing & Pringlell2007l) . and 
in both cases this suggests that fragmenting discs are marginally 
optically thick. Individual clumps (with higher densities) thus have 
optical depths much greater than unity, so the optically thick ap- 
proximation appears valid across a wide range of disc scales. It 
is clear, however, that the opacity is not strictly constant in real 
discs. Dust opacity is essentially independent of density, but varies 
strongly as a function of temperature: the opacity scales as for 
T < I60K, and is approximately constant for 160K< T < lOOOK 
( iBell & Lin|[T994h . It is reasonable to assume that the temperature 
will increase somewhat during collapse, leading to an increase in 
the dust opacity, so it seems likely that our adopted function D(p) 
somewhat over-estimates the efliiciency of radiative cooling. Thus 
cooling is likely slower than predicted by our simulations in real 
fragmenting clumps. Models of core-collapse during star forma- 
tion predict only moderate increases in temperature, however, so 
the uncertainly due to the temperature dependence of dust opac- 
ity is likely small. We note in passing that we also neglect further 
complexities of the collapse process, such as changes in gas com- 
position with increasing temperature (due to dust sublimation, H2 
dissociation, etc.): such issues are well beyond the scope of this 
work. 

An additional complication is that is still not clear whether or 



not it is possible to form a disc with /3 < 1 . A coherent disc structure 
cannot be assembled in the presence of such rapid cooling, as the 
growing disc would simply fragment on a dynamical time-scale. 
Various authors have argued in favour of rapid cooling in self- 
gravitating discs around massive black holes (e.g. lBertin & Lodatol 
[2001; Goodman 2003), but these arguments generally rely on the 
assumption of roughly steady-st ate cond i tions i n the disc, rather 
than considering disc formation. iRafikovl j2007h argued that con- 
vective cooling results in an cooling time of P - 0.1 , but again 
does not consid er disc form ation in detail. By contrast. iNavakshirj 
f2006, see also ILevirJl20a^ argues against the presence of rapid 
cooling in a fragmenting, star-forming disc around a massive black 
hole, instead arguing that the disc mass is accreted gradually. The 
cooling time therefore decreases slowly as the disc mass is accumu- 
lated, so fragmentation occurs at /3 ^ /Jcrit- If however, the disc was 
formed in a rapid "accretion event", it seems that fragmentation 
could occur with p < I . In addition the thermodynamics of the disc 
formation process are not fully understood, and one can conceive of 
a scenario where the disc is assembled with a longer cooling time 
and then begins to cool faster as it evolves (e.g. lciarke et al.l2007h . 
One possible mechanism appeals to the "opacity gap" at 10'-10'*K: 
if a disc was assembled from warm, > 1500K gas then the conden- 
sation of dust at ~ 1000-1500K would result in a dramatic increase 
in the cooling rate as the temperature dropped, potentially result- 
ing in a cooling time-scale shorter than the dynamical time-scale. 
Our simple model is not able to address this issue in detail, but this 
problem provides fertile ground for future research. 

Despite the fact that our simulations show good numerical 
convergence with increasing particle number, numerical resolution 
is still an issue that must be considered with care. In particular, 
as noted in Section [STI the dynamical interactions between differ- 
ent clumps are under-resolved even in our highest-resolution sim- 
ulations. This is a common feature of hydrodynamic simulations 
which seek to model large dynamic ranges in length and density, 
and in our simulations this results in unphysically rapid mergers of 
binary and multiple systems (as seen, for example, in FigO- We 
treat this problem by simply stopping our calculations before such 
under-resolved mergers begin to dominate the simulations, and at- 
tach only limited significance to our results beyond the point at 
which the first such mergers begin to happen. [Note, however, that 
real mergers between fragments are expected to occur in such sys- 
tems, and indeed have been invoked to explain the large discrep- 
ancy between the typical fragme nt mass and the observed stellar 
masses near the Galactic Center ( lLevinll2007h .] A more sophisti- 
cate d approach could , for example, make use of "s ink particles" 
(e.g.lBate etalj 19951) , modified thermal physics (e.g. Japps en et al.l 
I2OO5 ), or sub-resolution physics (e.g. Nayakshin et al. 200"^ in or- 
der to follow the calculations further forward in time, but in our 
models (and in particular in the simulations with fast cooling) it is 
not clear that such approaches would yield more accurate results. 
Here, a better approach may be to "re-simulate" a small region of 
our computational domain at much higher resolution in order to 
study the inter-clump dynamics in more detail, but this is beyond 
the scope of this paper. 



4.1 Application to the Galactic Centre 

We now consider the implications of our results for the forma- 
tion of stars near the Galactic Centre (GC). Many young, mas- 
sive stars are known to exis t close to the black hole at the cen- 
tre of the Galaxy (Ghez et al] |l998l : lGe"nzel et al.ll2003l : lGhez et all 
I2OO5I : IPaumard et al., aOQ^ . and many of the stars at =^ 0. Ipc 



© 2008 RAS, MNRAS QOCrTHTol 



Fragmenting discs with short cooling times 9 



are k nown to form a coherent ring or disc (e.g. iGenzel et al.l 
I2OO3IFI . A popular theory for the formation of these stars is 
fragmentatio n of an accretion disc due to gravitational insta- 
biUties (e.g. iLevin & Beloborodovl |2003| ; iGoodman & Taij|2004 
iNavakshin & CuadriJl2o65 !r and recent hydrodynamic simulations 
have found that the properties of the stars formed in this man- 
ner are broadly consis t ent with the observed G C stellar disc (e.g. 
INavakshin et alj l2007l : [Alexander et alj |2008|) . However, to date 
most of these simulations have used greatly simplified thermal 
physics, adopting scale-free cooling and/or simplified equations of 
state, and given the critical role of cooling in the development of 
gravitational instabilities this treatment is less than ideal. 

Additionally, for numerical reasons most previous simulations 
of the formation of the GC stellar disc have made use of relatively 
slow cooling time-scales (with =^ 3), although INavakshin et al.l 
1 I2OO7) did run a single model with = 0.3 (as discussed in Sec- 
tion|3]2}- However, analytic estimates of the cooling times in black 
hole discs generally find that the cooling time should be short com- 
pared to the dynamical time (e.g. Goodman 2003; Rafikov 2007) 
although, as discussed above, how such systems form remains un- 
certain. Estimates of the cooling parameter in a fragmenting disc 
at the GC range from /? ~ 0.1 l lRafikovllioOTh to ^6 = yScit ^ 3 
( lNavakshinl l2006^. which is approximately the range spanned by 
our simulations. As discussed above, we find that a dramatic de- 
crease in the cooling time leads to only a modest decrease in the 
characteristic stellar mass, and given the dependence of the charac- 
teristic mass on other disc properties (notably the disc aspect ratio 
H/R), this is unlikely to have a significant effect on the MF of the 
GC stellar disc. However, our simulations indicate that fast cooling 
results in a very high primordial binary/multiple fraction (of order 
unity), and although binaries do form when cooling is slower, they 
are less common. We therefore suggest that, if the GC stellar disc 
formed by fragmentation of a gaseous accretion disc, a high binary 
fraction may be indicative of rapid cooling in the parent accretion 
disc. 

It is not clear, however, if this prediction provides a means 
of distinguishing between different proposed scenarios for the for- 
mation of the GC stellar disc. The massive star population in the 
field is known to have a high multiplicity fraction, and the binary 
fraction for Wolf-Ra yet stars in the solar neighbourhood is ~ 40% 
( Ivan der Huchll200l[) . Little is currently known about binary statis- 
tics in the GC stella r cluster, with neith er observations nor dynam- 
ical arguments (e.g. ICuadra et al]|2008l) providing real constraints. 
It is entirely possible (and in fact likely) that many unresolved close 
binaries near the GC may be masquerading as single objects in cur- 
rent data, but some binary systems have been observed. In particu- 
lar IRS 16SW, an eclipsing binary with a period of ^19.5 days, ap- 
pears to be a contact binary consisting of two a pproximately equal- 
mass Wolf-Rayet s t ars with masse s > 50Mo dMartins et alj|2006l : 
IPeeples et al. I I2OO7I : iRafelski et al.1l2007l) . Systems of this type ap- 
pear to be a plausible outcome of fast cooling in a fragmenting, 
self-gravitating accretion disc, and we suggest that such systems 
should be common if the stellar disc at the GC formed via frag- 
mentation of a gaseous accretion disc. 



Note that there is some c ontroversy as to whe t her one or two discs of 
young stai-s exist at the GC JPaumard et al.ll2006t [Lu et alj 2006^. For the 
purposes of this discussion this is unimportant: the existence of one disc of 
young stars is sufficient to motivate our arguments, and they are not signifi- 
cantly altered by the existence of a second disc. 



5 SUMMARY 

We have used numerical hydrodynamic simulations to investigate 
the role of the cooling rate in the fragmentation of gravitationally 
unstable accretion discs. Our simulations make use of a simple 
cooling law, designed to mimic the behaviour of radiative cooling 
in an optically thick, geometrically thin disc. We performed de- 
tailed numerical convergence tests, and are satisfied that our simu- 
lations are well-resolved. We find that decreasing the cooling time 
results in the formation of a larger number of fragments, but with 
a correspondingly smaller characteristic mass. We also observed a 
high incidence of binaries, all of which form by dynamical cap- 
ture. Rapid cooling, on a time-scale shorter than the local dynam- 
ical time, results in a dramatic increase in the incidence of binary 
and multiple systems, although we are not able to determine the 
final properties of such systems. We have considered our results 
in the context of the young stellar disc(s) at the Galactic Centre, 
and find that - at the very least - observations of some stellar bi- 
naries close to the Galactic Center do not rule out formation in a 
fragmenting accretion disc. 
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